www.gusucode.com > Real Time Embedded Target Application of MPC 程序工具箱matlab源码 > Real Time Embedded Target Application of MPC/Chapter_9/Section_8/Chapter_9_Section_8_Script.m

    %Book Title: Practical Design and Application of Linear MPC
%Chapter: 9
%Section: 7
%Authors: Nassim Khaled and Bibin Pattel
%Last Modified: 10/07/2017
%%
clc
clear all
% Specify the path to look for the simulation results
curr_working_folder = pwd;
cd('..');
system_id_data_path = [pwd '\Section_6'];
cd(curr_working_folder);
% Load the system id data from Section 9.5
load_str = ['load ' system_id_data_path '\Chapter_9_Section_6_System_ID_Data.mat '];
eval(load_str);



%% Model # 1 Grouping simulation data and neglecting first few seconds(Nstart*stp_sz)
% Specify sample time of data
stp_sz = 0.2;
Nstart = 300;%50; % Clip the data starting from time t = 10 Secs to t = 145 Secs
Nend = 725;%
Measured_Outputs= Motor_Speed_RPM_Filtered(Nstart:Nend,2); % Measured outputs
Manipipulated_Variables= PWM_Input(Nstart:Nend,2); %Manipulated variables

% Obtaining initial conditions at step time=Nstart
Measured_Outputs_Nstart=Measured_Outputs(1,:); %Capturing the measured outputs at step time= Nstart
Manipipulated_Variables_Nstart=Manipipulated_Variables(1,:); %Capturing the manipulated variables at step time= Nstart

% Forcing response to start from zero initial condi-tions
Measured_Outputs_zero_initial_conditions=Measured_Outputs-repmat(Measured_Outputs_Nstart,length(Measured_Outputs),1); %Subtracting initial conditions for measured outputs to obtain zero response at step time= Nstart
Manipipulated_Variables_zero_initial_conditions=Manipipulated_Variables-repmat(Manipipulated_Variables_Nstart,length(Manipipulated_Variables),1); %Subtracting initial conditions for ma-nipulated variables to obtain zero actuation at step time= Nstart

% Prepare date for system identification
data=iddata(Measured_Outputs_zero_initial_conditions,Manipipulated_Variables_zero_initial_conditions,stp_sz); %data is packaged for system identification using idda-ta

% Generate a preliminary 2nd order system that fits the data
sys1=n4sid(data,2,'Form','canonical','DisturbanceModel','None','InputDelay',[0 0]','InitialState','zero'); %n4sid generates a preliminary system in the canonical form
%with zero disturbance, zero delay and zero initial conditions

% Generate a more refined system
sys2 = pem(data,sys1,'InitialState','zero') %pem gener-ates a more refined system that fits the data better

% Define the options for comparing the various identi-fied systems
opt = compareOptions('InitialCondition','Z');
[Y,fit,x0]=compare(data,sys2);
Y_1=Y.OutputData(:,1);

figure
subplot(2,1,1)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Manipipulated_Variables_zero_initial_conditions(:,1),'linewidth',2)
grid on
legend('PWM Command Value')
ylabel('PWM Command Value')
title('Motor Speed System Identification for Linear Region #1 Midpoint Step');
subplot(2,1,2)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Measured_Outputs_zero_initial_conditions(:,1),'linewidth',2)
grid on
hold on
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Y_1,'r','linewidth',2)
legend('Motor Speed Measured RPM','Motor Speed Modeled RPM')
ylabel('Motor Speed RPM')
xlabel('Time Seconds')
set(gcf,'color',[1 1 1]);

eval_str = 'Motor_Linear_Model_Mode_1 = sys2;';
eval(eval_str);
eval_str = ['save Motor_Linear_Model_Mode_1.mat' ' Motor_Linear_Model_Mode_1'];
eval(eval_str);
%% Model # 2 Grouping simulation data and neglecting first few seconds(Nstart*stp_sz)
Nstart = 1050;%800; % Clip the data starting from time t = 160 Secs to t = 245 Secs
Nend = 1475;%
Measured_Outputs= Motor_Speed_RPM_Filtered(Nstart:Nend,2); % Measured outputs
Manipipulated_Variables= PWM_Input(Nstart:Nend,2); %Manipulated variables

% Obtaining initial conditions at step time=Nstart
Measured_Outputs_Nstart=Measured_Outputs(1,:); %Capturing the measured outputs at step time= Nstart
Manipipulated_Variables_Nstart=Manipipulated_Variables(1,:); %Capturing the manipulated variables at step time= Nstart

% Forcing response to start from zero initial condi-tions
Measured_Outputs_zero_initial_conditions=Measured_Outputs-repmat(Measured_Outputs_Nstart,length(Measured_Outputs),1); %Subtracting initial conditions for measured outputs to obtain zero response at step time= Nstart
Manipipulated_Variables_zero_initial_conditions=Manipipulated_Variables-repmat(Manipipulated_Variables_Nstart,length(Manipipulated_Variables),1); %Subtracting initial conditions for ma-nipulated variables to obtain zero actuation at step time= Nstart

% Prepare date for system identification
data=iddata(Measured_Outputs_zero_initial_conditions,Manipipulated_Variables_zero_initial_conditions,stp_sz); %data is packaged for system identification using idda-ta

% Generate a preliminary 2nd order system that fits the data
sys1=n4sid(data,2,'Form','canonical','DisturbanceModel','None','InputDelay',[0 0]','InitialState','zero'); %n4sid generates a preliminary system in the canonical form
%with zero disturbance, zero delay and zero initial conditions

% Generate a more refined system
sys2 = pem(data,sys1,'InitialState','zero') %pem gener-ates a more refined system that fits the data better

% Define the options for comparing the various identi-fied systems
opt = compareOptions('InitialCondition','Z');
[Y,fit,x0]=compare(data,sys2);
Y_1=Y.OutputData(:,1);

% Plot the Model Fit
figure
subplot(2,1,1)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Manipipulated_Variables_zero_initial_conditions(:,1),'linewidth',2)
grid on
legend('PWM Command Value')
ylabel('PWM Command Value')
title('Motor Speed System Identification for Linear Region #2 Midpoint Step');
subplot(2,1,2)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Measured_Outputs_zero_initial_conditions(:,1),'linewidth',2)
grid on
hold on
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Y_1,'r','linewidth',2)
legend('Motor Speed Measured RPM','Motor Speed Modeled RPM')
ylabel('Motor Speed RPM')
xlabel('Time Seconds')
set(gcf,'color',[1 1 1]);

eval_str = 'Motor_Linear_Model_Mode_2 = sys2;';
eval(eval_str);
eval_str = ['save Motor_Linear_Model_Mode_2.mat' ' Motor_Linear_Model_Mode_2'];
eval(eval_str);
%% Model # 3 Grouping simulation data and neglecting first few seconds(Nstart*stp_sz)
Nstart = 1800;%1550; % Clip the data starting from time t = 160 Secs to t = 245 Secs
Nend = 2225;%2475;%
Measured_Outputs= Motor_Speed_RPM_Filtered(Nstart:Nend,2); % Measured outputs
Manipipulated_Variables= PWM_Input(Nstart:Nend,2); %Manipulated variables

% Obtaining initial conditions at step time=Nstart
Measured_Outputs_Nstart=Measured_Outputs(1,:); %Capturing the measured outputs at step time= Nstart
Manipipulated_Variables_Nstart=Manipipulated_Variables(1,:); %Capturing the manipulated variables at step time= Nstart

% Forcing response to start from zero initial condi-tions
Measured_Outputs_zero_initial_conditions=Measured_Outputs-repmat(Measured_Outputs_Nstart,length(Measured_Outputs),1); %Subtracting initial conditions for measured outputs to obtain zero response at step time= Nstart
Manipipulated_Variables_zero_initial_conditions=Manipipulated_Variables-repmat(Manipipulated_Variables_Nstart,length(Manipipulated_Variables),1); %Subtracting initial conditions for ma-nipulated variables to obtain zero actuation at step time= Nstart

% Prepare date for system identification
data=iddata(Measured_Outputs_zero_initial_conditions,Manipipulated_Variables_zero_initial_conditions,stp_sz); %data is packaged for system identification using idda-ta

% Generate a preliminary 2nd order system that fits the data
sys1=n4sid(data,2,'Form','canonical','DisturbanceModel','None','InputDelay',[0 0]','InitialState','zero'); %n4sid generates a preliminary system in the canonical form
%with zero disturbance, zero delay and zero initial conditions

% Generate a more refined system
sys2 = pem(data,sys1,'InitialState','zero') %pem gener-ates a more refined system that fits the data better

% Define the options for comparing the various identi-fied systems
opt = compareOptions('InitialCondition','Z');
[Y,fit,x0]=compare(data,sys2);
Y_1=Y.OutputData(:,1);

% Plot the Model Fit
figure
subplot(2,1,1)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Manipipulated_Variables_zero_initial_conditions(:,1),'linewidth',2)
grid on
legend('PWM Command Value')
ylabel('PWM Command Value')
title('Motor Speed System Identification for Linear Region #3 Midpoint Step');
subplot(2,1,2)
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Measured_Outputs_zero_initial_conditions(:,1),'linewidth',2)
grid on
hold on
plot(0:stp_sz:stp_sz*(length(Y_1)-1),Y_1,'r','linewidth',2)
legend('Motor Speed Measured RPM','Motor Speed Modeled RPM')
ylabel('Motor Speed RPM')
xlabel('Time Seconds')
set(gcf,'color',[1 1 1]);

eval_str = 'Motor_Linear_Model_Mode_3 = sys2;';
eval(eval_str);
eval_str = ['save Motor_Linear_Model_Mode_3.mat' ' Motor_Linear_Model_Mode_3'];
eval(eval_str);